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PREFACE 


This  is  one  of  a series  of  reports  on  research  and  development 
being  carried  out  under  the  Rand/ARPA  Climate  Dynamics  Project.  While 
other  reports  have  given  research  or  experimental  results,  the  present 
report  is  aimed  primarily  at  development,  specifically  model  improve- 
ment in  one  important  area:  the  need  for  an  interactive  world  ocean 

model  for  climate  experiments  in  conjunction  with  the  Mintz-Arakawa 
atmospheric  general  circulation  model.  This  report  concentrates  on 
formulation  of  the  basic  ocean  model,  with  application  to  the  Pacific. 
Generalization  to  the  world  ocean  and  coupling  with  the  atmospheric 
counterpart  will  be  treated  in  separate  publications. 

At  least  two  previous  Rand  publications  provided  direct  inputs  to 
the  present  work: 

RM-6211-ARPA,  Numerical  Studies  of  Planetary  Circulations  in  a 
Wind-Driven  Ocean  on  the  Sphere. 

X-505-ARPA,  A Calibrated  Analytical  Model  for  the  Thermohaline 
and  Wind-Driven  Circulation  in  the  Interior  of  a Subtropical 
Ocean. 

Other  Rand  publications  on  ocean  circulations  include: 

RM-5594-XSF , The  Ekman  Vertical  Velocity  in  an  Enclosed  S-Plane 
Ocean. 

RM-6110-RC,  The  Effects  of  Western  Coastal  Orientation  on  Rossby- 
Wave  Ref  ection  and  the  Resulting  Large  Scale  Oceanic  Circulation. 

RM-62 10-ARPA,  A Note  on  the  Lateral  Eddy  Viscosity  Due  to  Transient 
Rossby  Waves  in  a Barotropic  Model. 

The  Climate  Dynamics  Project  is  sponsored  bv  the  Defense  Advanced 
Research  Projects  Agency. 
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SUMMARY 


PURPOSE 

Ihis  research  was  undertaken  to  develop  a numerical  model  for 
monthly  and  seasonal  predictions  of  temperatures  and  velocities  in  the 
upper  layers  of  the  ocean.  The  Pacific  model  represents  one  phase  in 
the  development  of  a world  ocean  model  for  climate  dynamics  experiments 
with  a coupled  atmospheric  model. 

numerical  MODEL 

The  circulation  and  temperature  distribution  in  the  upper  layers 
ot  the  Pacific  Ocean  are  simulated  bv  a two-level  numerical  model  with 
a horizontal  bottom  at  300  m depth.  The  model  is  based  on  the  primitive 
Eulerian  equations  and  contains  a free  upper  surface.  The  motion  is 
driven  by  prescribed  distributions  of  surface  heating  and  wind  stress, 
and  is  retarded  by  frictional  stresses  at  the  bottom  and  at  the  coasts. 
An  explicit  numerical  scheme  is  used  which  employs  central  space-  and 
time-differences  with  the  exception  that  the  diffusion  terms  are  lagged 
one  time  step  for  stability. 

Solutions  for  temperature  and  velocity  are  sought  as  an  initial 
value  problem.  Initial  conditions  are  nondivergent  and  gcostrophic  for 
prescribed  initial  temperature  fields.  Initialization  is  completed  by 
barotropic  spinup  of  the  model  while  holding  the  initial  temperature 
and  wind  stress  fixed.  Predictions  can  then  be  made  by  introducing  the 
temperature  prediction  equation  together  with  surface  heating,  and  by 
allowing  temporal  variations  of  wind  stress. 

DEMONSTRATION  RUN 

The  demonstration  run  reported  here  employs  January  climatological 
distributions  of  temperature,  wind  stress,  and  heating  as  initial  and 
boundary  data.  Initialization  gives  strong  flow  in  the  equatorial  re- 
gion having  a velocity  reversal  with  depth,  in  rough  agreement  with 
observation.  However,  the  north-south  components  are  too  large  because 
of  an  inaccurate  specification  of  initial  temperature  in  low  latitudes, 
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where  data  are  sparse.  A 15-day  prediction  with  the  model  shows  the 
simultaneous  development  of  more  realistic  temperatures  and  velocities 
in  the  equatorial  zone,  including  a relatively  cold  geographic  equator 
with  colder  water  toward  the  east,  and  verv  nearly  zonal  flow  directed 
westward  at  the  upper  level,  eastward  at  the  lower  level.  However, 
countercurrents  displaced  north  or  south  of  the  equator  are  not  dis- 
cerned in  the  results,  and  some  features  of  the  flow  near  eastern 
boundaries  may  be  unrealistic. 

Although  the  demonstration  run  reported  here  represents  only  a 
special  case,  the  results  found  in  the  equatorial  zone  are  typical  of 
results  found  from  many  experiments  made  with  this  and  simpler  versions 

of  the  model. 

CONCLUSIONS 

• A main  conclusion  from  the  point  of  view  of  climate  dynamics 
is  that  an  oceanic  general  circi’1  ution  model  can  make  15-day 
predictions  showing  appreciable  changes  in  temperature  dis- 
tributions, at  least  in  low  latitudes.  This  may  have  im- 
portant implications  for  monthly  and  seasonal  simulations 
with  interacting  atmospheric  and  oceanic  models. 

• Deficiencies  in  the  results  can  for  the  most  part  be  traced 
to  known  limitations  of  the  model  design.  Tnese  include 
problems  of  resolution  near  boundaries  and  across  the  equa- 
tor, and  the  lack  of  continental  shelf  topography.  The 
limitations  are  probably  not  serious  in  light  of  the  limita 
tions  of  existing  atmospheric  general  circulation  models. 

• The  results  are  generally  encouraging,  largely  because  of 
the  realism  of  the  predicted  large-scale  temperature  anomaly 
in  the  eastern  tropics.  Moreover,  computer  times  are  modest 
compared  with  those  required  for  a two-level  atmospheric 
model.  The  Pacific  model  is  being  extended  to  the  world 
ocean  case,  with  possible  improvements  postponed  to  later. 

• From  an  oceanographic  viewpoint,  the  results  indicate  that 
the  occurrence  of  zonal  flow  near  the  equator  '1s  intimately 


related  to  the  thermal  structure.  In  particular,  an  east- 
west  temperature  gradient  is  required  by  the  model.  However, 
a physical  explanation  is  not  offered  here,  and  further  study 
with  a more  specialized  model  is  required  to  clarify  this 
point. 
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SYMBOLS 


( ) = 

*h  = 

A = 
ra 

H = 
K - 


a vector  quantity  ( ) is  underscored  with  a tilde 
constant  horizontal  eddy  thermal  conductivity 
constant  horizontal  eddy  kinematic  viscosity 
constant  total  mean  depth 
bottom  friction  coefficient 

total  vertical  heat  flux  at  the  air-sea  interface,  positive 
downward 

total  vertical  heat  transport  at  level  k 
temperature 

temperature  at  level  k 

horizontal  velocity  vector,  vertically  averaged  for  layer  k 
eastward,  northward  components  of  U 


a = earth's  mean  radius 

cp  = specific  heat  of  sea  water  at  constant  pressure 
f = coriolis  parameter,  2Q  sin  cp  b (u  tan  cp)  / a 
8 = gravitational  acceleration 
hk  = constant  mean  thickness  of  layer  k 

i,j,k  = discrete  indices  denoting  eastward,  northward,  and  downward 
position 

- ~ un^t  vertical  vector,  positive  upward 
n - unit  normal  vector,  positive  outward 
P = pressure 

Pa  = atmospheric  pressure 

q - penetrating  solar  radiation  source  term  in  heat  flux 
t = time 

U.v.w  = eastward,  northward,  and  upward  velocity  components 
v = horizontal  velocity  vector 

Wk  = vertical  velocity  at  level  k,  positive  upwards 
Ws  = vertical  velocity  at  the  air-sea  interface 

z = vertical  coordinate,  positive  upward  with  respect  to  mean 
sea  level 


x = coefficient  of  thermal  expansion 
At  = time  increment 

, <P  = east-west  and  north-south  grid  spacing 

elevation  of  air-sea  interface  with  respect  to  mean  sea  level 
k = vertically  averaged  temperature  for  layer  k 
< " constant  vertical  eddy  thermal  conductivity 
1 = longitude 

i discrete  longitudinal  location  of  grid  points 
v = constant  vertical  eddy  kinematic  viscosity 
\ = vertically  averaged  pressure  for  layer  k 
P = density  at  constant  pressure  and  salinity 
= reference  density 
Ifc  ~ bottom  stress  vector 
Is  = winc*  stress  vector 

•k 

Tk  « total  vertical  transport  of  horizontal  momentum  at  level  k 

Cp  _ 

- eastward,  northward  components  of  stress  vectors 
<p  = latitude 

cpj  * discrete  latitudinal  location  of  grid  points 
~ earth's  rotation  rate 


I . INTRODUCTION 


The  world  ocean  forms  a vast  and  important  connecting  link  to  the 
atmospheric  heat  engine,  and  hence  to  climatic  changes,  via  exchanges 
of  heat,  momentum,  and  water  vapor  across  the  air-sea  interface.  Among 
the  many  interactions  on  many  scales  of  motion  and  time,  that  of  "wind- 
stress  driving/sea-temperature  changing"  appears  crucial  for  more  ac- 
curate prediction  of  large-scale  circulations  resulting  from,  say,  some 
environmental  perturbation.  Changes  in  atmospheric  wind  stresses  acting 
on  the  ocean  cause  changes  in  transports  of  heat  in  the  latter,  which 
changes  the  sea  temperature.  This  changes  the  distributions  of  heat  and 
water  vapor  fluxes  between  ocean  and  atmosphere  which,  in  turn,  changes 
atmospheric  wind  patterns.  Although  this  is  understood  in  a rough, 
qualitative  way,  more  precise  information  is  needed  in  studies  of  cli- 
matic changes.  The  cause  of  the  changed  wind  stress  pattern  may  have 
its  origin  in  some  remote  region  of  the  globe;  we  need  to  know  how  much 
the  sea  temperature  changes,  and  where  and  how  soon,  and  we  need  to 
know  more  precisely  the  effect,  in  turn,  on  the  atmosphere. 

One  approach  to  studying  such  large-scale  climatic  effects  and 
interactions  is  by  numerical  modeling  of  the  combined  atmosphere-ocean 
system.  Numerical  models  of  the  global  atmosphere  are  already  well 
established  and  operational.  The  better  known  models  in  this  country 
include  those  of  the  Geophysical  Fluid  Dynamics  Laboratory  (Smagorinsky 
et  al.  1965),  the  National  Center  for  Atmospheric  Research  (Kasahara 
and  Washington  1971),  and  the  University  of  California  at  Los  Angeles 
(Mintz  1968,  Gates  et  al.  1971).  Sea  temperatures  in  these  models  have 
been  specified  climatologically , however,  and  there  has  been  no  fully 
interacting  ocean. 

The  development  of  comparable  ocean  models  began  somewhat  later, 
notably  with  the  work  of  Bryan  and  Cox  (1967,  1968),  and  has  not  yet 
matured.  The  process  of  coupling  an  oceanic  model  to  its  atmospheric 
counterpart  presents  further  problems.  Pioneering  efforts  in  this 
direction  have  been  made  by  Manabe  and  Bryan  (Manabe  1969),  who  studied 
a simple  interacting  model  of  atmospheric  flow  over  a single  rectangular 
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ocean  and  a comparably  idealized  land  distribution.  More  recently, 

Bryan  and  Cox  (1972)  have  presented  results  from  their  noninteracting 
world  ocean  model  for  the  case  of  a homogeneous  ocean. 

In  their  pioneering  ocean  modeling,  Bryan  and  Cox  placed  more 
emphasis  on  studying  the  ocean  in  its  own  right  (and  properly  so)  over 
long  periods  of  ocean  time,  and  perhaps  less  on  developing  a somewhat 
cruder  ocean  model  for  routine  use  in  studying  the  joint  atmosphere- 
ocean  system.  Changes  occur  very  slowly  in  the  deep  ocean,  with  time 
periods  measurable  in  centuries.  To  make  such  time-integrations  feas- 
ible on  existing  computing  machines,  Bryan  and  Cox  adopted  a rigid  lid 
approximation  in  which  the  vertical  velocity  is  set  identically  zero 
at  the  (level)  sea  surface.  This  effectively  removes  surface  gravity 
waves,  apparently  without  appreciably  affecting  the  resulting  steady 
circulation  patterns,  and  thereby  allows  time  increments  to  be  orders 
of  magnitude  larger  than  would  otherwise  have  been  possible  with  simple 
explicit  numerical  schemes. 

This  method  is  not  without  a price  in  terms  of  computational  com- 
plexity. Pressure  can  no  longer  be  evaluated  directly  from  the  hydro- 
static equation,  because  of  an  unknown  reference  pressure  (e.g. , pressure 
at  the  rigid  lid  can  no  longer  be  specified  as  atmospheric  pressure), 
and  is  usually  eliminated  from  the  horizontal  momentum  equations  by 
cros s— di f f erent iation . This  raises  the  mathematical  order  of  the  equa- 
tions and  gives  a Poisson  boundary  value  problem  to  be  solved  at  every 
ime  step  for  a stream  function.  The  (mathematical)  connectivity  of 
the  solution  region  becomes  an  important  question  in  such  problems  for 
reasons  of  uniqueness.  Introduction  of  a single  island  into  a simple 
closed  basin  causes  the  solution  region  to  become  doubly  connected,  and 
makes  the  problem  more  difficult  to  solve.  The  multiple  connectivity 
of  the  world  ocean  adds  to  difficulties  of  this  sort.  Nevertheless, 
long-term  integrations  have  proved  feasible  by  this  method,  at  least 
for  simple  closed  basins  (Bryan  and  Cox  1968,  Bryan  1969b,  Cox  1970), 
and  much  has  been  learned  from  these  solutions. 

However,  a slight  inconsistency  may  arise  when  an  ocean  model 
filtered  with  the  rigid  lid  approximation  is  coupled  to  an  atmospheric 
model  based  on  the  unfiltered,  primitive  Eulerian  equations  (that  is, 
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unf iltered  in  terms  of  external  gravity  waves).  Most  of  the  known 
atmospheric  models  are  apparently  of  the  unfiltered  type.  For  these 
models,  integrations  spanning  centuries  would  not  be  feasible  with 
existing  computers.  Bryan  (1969b)  and  Manabe  (1969)  have  formulated 
a promising  approach  which  may  lead  to  fruitful  investigations  of 
long-term  border  of  centuries)  changes  in  global  climate.  The  essence 
of  the  method  is  that  the  atmospheric  model  is  integrated  for  a rela- 
tively short  period  uniformly  incremented  over  the  course  of  the  long- 
term integration  of  the  oceanic  ■ ^nterpart.  Some  justificati  a for 
this  is  found  in  the  relatively  rapid  response  of  the  atmosphere  com- 
pared with  that  of  tht  ocean. 

Conversely,  if  jointly-running  atmospheric  and  oceanic  models  are 
to  simulate  the  same  time  span,  then  such  times  are  apparently  limited 
to  several  years,  a decade  perhaps,  with  present  technology.  This  is 
in  an  equally  important  area  of  seasonal  changes  in  climate.  Of  course 
a f iltered  ocean  model  can  (and  no  aoubt  will)  be  used  in  such  investi- 
gations. However,  the  need  f - larger  time  increments  becomes  less 
pressing,  and  the  disadvantage  of  the  filtered  models  may  make  the 
simplicity  of  the  unfiltered  models  more  appealing. 

Thus,  there  appears  to  be  a distinct  need  for  an  ocean  model  which 
may  be  simpler  in  concept  and  design,  and  which  is  intended  only  for 
simulations  of  from  a few  weeks  duration  to  a few  year's. 

This,  then,  is  a resumption  of  the  work  begun  by  Crowley  (1968) 
to  construct  a model  of  the  stratified  circulation  of  the  world  ocean 
based  on  the  unfiltered  primitive  equations.  Crowley's  model  may  have 
been  several  years  ahead  of  its  time;  a new  generation  of  computers 
that  allow  minimal  (spatially  uniform)  grid  resolution  of  Gulf-Stream- 
type  phenomena  in  the  world  ocean  is  only  now  becoming  a reality.  The 
present  formulation  will  also  take  advantage  of  some  simple  concepts 
of  conservation  of  quadratic  quantities  for  nonlinear  stability  (al- 
though not  in  an  exact  sense  in  the  piesent  model),  as  discussed  by 
Bryan  (1966). 

In  many  respects,  with  one  notable  exception,  the  formulation  is 
a generalization  of  the  homogeneous  ocean  model  of  Gates  (1968).  The 
exception  is  that,  unlike  Gates'  model,  the  present  model  neglects 
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bottom  topography.  The  purpose  here  is  to  approximate  the  thermally 
active  circulation  in  only  the  upper  few  hundred  meters  of  the  ocean 
for  seasonal  simulations.  The  model,  again  like  Crowley's,  contains 
a horizontal  bottom  and  vertical  coastlines.  But  unlike  Crowley's 
model,  which  contained  six  levels  extending  to  2000  m depth,  the  present 
model  contains  only  two  levels,  and  the  total  depth  is  300  m.  Inspira- 
tion, and  some  justification,  for  a two-level  ocean  model  is  found  in 
the  relative  success  achieved  with  simple  two-level  atmospheric  models. 
Indeed,  the  plan,  and  hope,  is  to  develop  an  ocean  model  that  might  be 
a true  companion  to  the  Rand  version  of  the  Mintz-Arakawa  two-level 
atmospheric  circulation  model. 

This  report  concentrates  on  formulation  of  the  basic  model  for  a 
bounded  ocean  basin,  with  application  to  the  Pacific.  Generalization 
to  the  world  ocean,  and  coupling  with  the  atmospheric  counterpart,  will 
be  treated  in  separate  publications. 


II.  PHYSICAL  AND  MATHEMATICAL  MODEL 


A.  PHYSICAL  MODEL 

We  wish  to  model  the  upper  layers  of  the  ocean  in  order  to  consider 
seasonal  effects  in  the  general  circulation.  Longer-term  effects  asso- 
ciated with  the  deeper  circulation  are  neglected.  Therefore,  consider 
as  a model  a thermally  stratified  fluid  which  fills  a shallow  basin  of 
uniform  depth  H on  a rotating  spherical  earth.  The  fluid  is  bounded 
laterally  by  continental  coastlines,  assumed  vertical.  Coasts  of 
islands  likewise  are  assumed  vertical.  For  definiteness,  the  depth  H 
will  be  taken  as  300  m,  which  corresponds  roughly  to  the  observed  base 
of  strong  subsurface  currents  in  equatorial  regions  (Knauss  1960). 

In  the  simplest  case,  which  is  treated  here,  the  fluid  is  thermally 

stratified,  and  salinity  stratification  is  neglected.  Again  for  defin- 
iteness, the  flud  is  taken  to  be  of  uniform  constant  sn'inity,  say  35°/ 

(parts  per  thousand)  for  average  sea  water,  and  obeys  a simple  linear 

equation  of  state  relating  the  density  at  atmospheric  pressure  to  the 
temperature.  This  equation  will  be  a linearization  based  on  the  formula 
proposed  by  Eckart  (1958).  Changes  of  phase  will  not  be  allowed  in  the 
model.  Any  sea  ice  present  in  high  latitudes  will  be  specified  (from 
observation)  in  the  model  by  arbitrarily  setting  the  sea  surface  temper- 
ature equal  to  the  (laboratory)  freezing  point  of  saline  water  of  35°/ 
(~-1.9°C).  °° 

ine  motion  is  driven  by  fluxes  of  heat  and  momentum  across  a free 
upper  surface,  and  is  retarded  by  turbulent  frictional  stresses  at  the 
coasts  and  possibly  also  at  the  bottom.  All  solid  boundaries  are 
assumed  insulated  and  impermeable.  The  upper  surface  heat  and  momentum 
(wind  stress)  fluxes  are  specified  climatologically  as  functions  of 
position  and  of  time.  In  general,  the  surface  heat  flux  may  inc1  jt 
contributions  proportional  to  the  air-sea  temperature  difference,  in 
which  case  the  surface  air  temperature  will  be  specified  climatologically. 
Haney  (1971b)  has  presented  a simple  and  lucid  formulation  of  a surface 
heat  flux  boundary  condition  for  ocean  circulation  models.  Formulation 
of  the  model  directly  in  terms  of  heat  fluxes  as  well  as  wind  stresses 
may  facilitate  coupling  the  ocean  model  to  its  atmospheric  counterpart. 
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One  feature  of  the  present  model  which  may  be  open  to  serious 
question  is  the  assumption  of  a horizontal  bottom  and  abrupt,  precipitous 
coastlines,  i.e.,  the  neglect  of  bottom  topography.  Evidence  presented 
by  Sarkisyan  and  Ivanov  (1971)  and  others  suggests  that  total  (vertically 
integrated)  transports  can  be  significantly  different  in  a stratified 
(or  baroclinic)  ocean  with  realistic  bottom  topography  as  opposed  to 
a stratified  ocean  with  a horizontal  bottom.  However,  it  is  less  clear 
how  topography  may  affect  upper  layer  transports  in  a stratified  ocean. 
Nevertheless,  one  of  several  possible  future  improvements  envisaged  for 
the  present  model  is  to  parameterize  a nonvanishing  vertical  vt  rocity 
at  the  lowermost  level  in  the  shalluw  ocean  model,  to  partially  account 
for  topographic  effects.  In  the  interim,  investigation  of  the  simpler, 
horizontal-bottom  model  seems  worthwhile. 

Section  VII  further  discusses  the  limitations  of  the  model. 

B.  GOVERNING  DIFFERENTIAL  EQUATIONS 

Large-scale  ocean  circulations  are  governed  by  the  primitive 
Eulerian  equations  for  an  incompressible,  inhomogeneous  fluid.  The 
momentum  equations  are  in  nearly  complete  form  in  the  horizontal,  but 
hydrostatic  in  the  vertical,  and  the  Boussinesq  approximation  is  em- 
ployed. The  equations  are  written  in  a spherical  coordinate  system 
with  \ and  c p denoting  longitude  and  latitude  respectively.  The  vertical 
coordinate  is  z,  taken  positive  upward  from  mean  sea  level.  Undiffer- 
entiated factors  (a  + z)  = (radius)  in  the  flux-divergence  terms 
are  approximated  by  a , where  a is  earth's  mean  radius.  In  vector 
notation  and  flux-divergence  form,  the  equations  are 

3v  a 3v 
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p ■ Poo  (1  - aT)  (2.4) 

It  = v • ^ - yt>  + It  « If  - wT  + q)  (2.5) 

where  v = (u,  v)  is  the  horizontal  velocity  vector  and  V the  horizontal 
differential  operator.  The  term  V • vv  in  Eq . (2.1)  denotes  a vector 
having  components  7 • (vu)  and  V • (w) . The  gradient  and  divergence 
terms  are  defined  by 
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The  horizontal  eddy  coefficients  of  momentum,  A , and  heat,  A,  , will 

4 m _ n 

be  assumed  constant.  Thus,  for  instance,  V • A Vv  = A V v. 

m ~ m 

In  Eq.  (2.5)  the  heat  flux  term  q is  given  by  S/p  c , where  S is 

OC  p 

short-wave  radiation  penetrating  downward  across  the  air-sea  interface. 
This  term  will  not  be  treated  explicitly  in  the  following;  its  contri- 
bution will  be  absorbed  into  the  total  surface  heat  flux  Qg , and  vertical 
integrations  will  be  taken  over  depths  greater  than  those  (~  50  m)  to 
which  appreciable  amounts  of  solar  radiation  can  penetrate. 

The  horizontal  bottom  of  the  model  ocean  is  placed  at  z = - H,  ard 
the  upper  surface  is  at  z = ?(A,  cp,  t).  It  is  assumed  that  |c|  « H. 

A vertical  section  of  the  model  near  a coastline  is  shown  schematically 
in  Fig.  1,  in  which  amplitudes  of  £ are  grossly  exaggerated,  and  in 
which  certain  other  quantities  appea  that  will  be  defined  presently. 

Boundary  conditions  for  the  model  are  as  follows.  At  the  air-sea 
interface,  z = C (X , cp,  t): 


bounded 
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where  x is  the  horizontal  wind  stress  vector  and  Q is  the  total  ver- 
tical  heat  flux  across  the  air-sea  interface.  Equations  (2.7)  and  (2.8) 
represent  a slight  approximation  consistent  with  the  hydrostatic  form 
of  the  governing  equations.  In  a more  general  formulation,  Eq . (2.7), 
for  instance , would  equate  the  normal  derivative  of  the  tangential  (as 
opposed  to  horizontal)  component  of  velocity  evaluated  at  z - ^ to  the 
tangential  component  of  wind  stress.  This  expression  hardly  differs 
from  that  given  above  if  sea  surface  slopes  are  small. 

At  the  bottom,  z = - H: 


w = 0 

(2.9) 

3y  i 

lb 

(2.10) 

3T 

^ = q = ° 

(2.11) 

where  an  assumption  is  required  to  relate  the  bottom  stress  T.  to  the 

~b 

velocity  field.  The  simplest  assumption  consistent  with  Fig.  1 is  that 
Tb  = 0.  A slightly  more  general  formulation,  which  is  used  in  the  model, 
assumes  that  the  bottom  stress  is  proportional  to  the  average  velocity 
in  the  lower  layer.  The  "friction"  coefficient  of  proportionality  may 
then  be  set  equal  to  zero  to  give  the  simpler,  no-stress  case. 

Conditions  at  the  coasts  are 

v = (normal  component  of  VT)  = 0 (2.12) 

Integration  of  the  continuity  Eq.  (2.3)  between  the  bottom  and  the 
sea  surface,  and  application  of  the  boundary  conditions  (2.6)  and  (2.9) 
gives  the  first  of  the  following  two  equations: 
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Leibnitz's  rule  was  employed  to  combine  the  integrated  divergence  with 
the  advection  term  in  arriving  at  (2.13).  The  pressure  Eq . (2.14)  is 
obtained  by  combining  (2.2)  and  (2.4)  and  integrating  between  a (vari- 
able) level  z and  the  sea  surface  where  p is  atmospheric  pressure. 

There  appears  to  be  no  exact  energy  integral  for  the  foregoing 
system  of  governing  equations  when  the  fluid  has  a free  surface.  This 
difficulty  is  not  encountered  in  models  in  which  ail  boundaries  are 
rigid,  because  the  processes  of  time-differentiation  and  integration 
over  a fixed  volume  can  be  freely  interchanged.  The  difficulty  arises 
here  principally  because  of  an  ambiguity  in  the  definition  of  density, 

0 versus  o,a.  Stated  another  way,  the  Boussinesq  approximation  should 
be  applied  after  the  more  exact  expression  for  energy  is  derived.  Thus, 
if  0^  is  replaced  with  0 in  Eq . (2.1)  and  energy  is  derived  in  the 
usual  way  except  that  integration  is  over  the  entire  nass  of  fluid  (as 
opposed  to  volume  integration),  then  the  order  of  time-differentiation 
and  integration  can  be  interchanged.  There  results: 


dt  (K  + P) 
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Here,  the  mass  integrals  have  been  changed  back  to  volume  integrals  by 
Cue  transformation  d(mass)  = .dxdydz.  W represents  the  volume  integral 
of  the  work  dene  by  the  eddy  diffusion  terms  in  (2.1),  Q the  heating 
done  by  the  eddy  conduction  terms  and  the  source  q in  (2.5).  The  sur- 
face integral  represents  the  work  done  by  atmospheric  pressure,  where 
v3  is  the  three-dimensional  velocity  vector  and  n the  u^t  outward  nor- 
mal to  tne  sea  surface.  The  above  has  the  usual  integral  K . + P = const 
in  the  absence  of  eddy  diffusion  processes,  sources  of  heat,  ar.d  work 
done  bv  atmospheric  pressure. 

Consistent  with  the  Boussinesq  approximation,  p is  now  approxi- 
mated by  in  the  kinet Lc  energy  (2.15)  but  is  retained  in  its  more 
complete  form  in  (2.16). 

To  this  point  no  additional  approximations  have  been  introduced 
beyond  those  implicit  in  the  hydrostatic  and  Boussinesq  form  or  the 
governing  equations.  The  foregoing,  except  for  Eq . (2.13),  is  similar 
to  the  system  of  equations  used  by  Crowley  (1968),  to  which  the  reader 
may  wish  to  refer  for  a derivation  from  the  more  complete  three- 
dimensional  equations  in  a spherical  coordinate  system.  Use  of  Eq. 
(2.13)  for  predicting  surface  elevation  was  pointed  out  by  Gates  (1968). 

The  foregoing  is  to  be  solved  numerically  as  an  initial  value 
problem.  For  given  initial  distributions  of  temperature,  velocity,  and 
surface  elevation,  predictions  of  these  same  quantities  are  to  be  made 
subject  to  the  boui  darv  conditions,  including  prescribed  sea-surface 
distributions  of  heating  and  wind  stress.  Evaluations  of  total  energy 
during  the  course  of  the  integration  provide  a measure  of  the  stability 
of  the  solution,  and  of  the  approach  to  a steadv  state. 
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III.  FINITE  DIFFERENCE  EQUATIONS 

. APPROXIMATED  LAYER  EQUATION’S : VERTICAL  DIFFERENCING 

Vertical  differencing  requires  special  care  for  several  reasons, 
and  is  therefore  treated  separately.  One  principal  reason  is  the 
hydrostatic  assumption  relating  to  the  vertical  coordinate.  Another 
is  the  relative  crudeness  of  resolution  of  a model  possessing  only  two 
levels  in  the  vertical.  Some  care  is  required  to  justify  what  may 
appear  to  be  rather  overly  detailed  vertical  distributions  of  the  de- 
pendent variables  in  a two-level  model.  The  approach  is  ty  means  of 
vertical  integration.  However,  the  results  will  be  formally  the  same 
as  differential/vertical-difference  equations. 

We  approximate,  first  by  . etaining  tee  surface  elevation  £ onlv 
to  lowest  order  in  Eqs.  (2.13)  and  (2.14).  The  upper  limits  of  inte- 
gration are  replaced  with  zeros,  and  the  resulting  system  of  equations 
is  treated  as  if  the  upper  surface  were  horizontal.  In  (2.13)  this 
amounts  to  ..eglecting  the  divergence  in  the  layer  of  thickness  lc|  com- 
pared with  the  divergence  in  the  much  larger  layer  of  constant  thick- 
ness H.  (The  sign  of  the  contribution  is  changed  for  £ < 0.)  The 
approximation  will  permit  interchanging  the  order  of  differentiation 
and  integration,  and  ail  resulting  variables  except  the  nonlinear  ad- 
vecticns  can  be  written  directly  in  terms  of  their  integrated  averages. 

Secondly,  in  the  horizontal  advection  terms  we  will  approximate 
the  average  of  products  with  the  product  of  averages,  e.g.,  vT  trvT. 

This  approximation,  although  common  in  numerical  schemes,  is  more  dif- 
ficult to  justify,  especially  in  a crude  two-level  model.  Some  justifi- 
cation is  found  in  the  quadrat ic-consc l v 'ng  (and  hence  stabilizing) 
properties  of  the  numerical  scheme  when  simple  averages  of  variables 
are  used.  The  approach  used  here  in  treating  the  advection  terms  is 
described  bv  Bryan  (1966). 

Define  two  layers  of  constant  thickness  h^  (layer  1)  and  h ^ (laver 
3)  as  shown  in  Fig.  1.  Layer  1 is  bounded  by  level  0 (mean  sea  level) 
and  the  intermediate  level  2;  layer  3 is  between  l">vel  2 and  level  4 
at  the  horizontal  bottom.  Vertically  averaged  values  of  v and  T for 
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layer  k (k  = 1,  3)  are  denoted  by  and  0^  as  shown  in  Fig.  1,  and 
are  defined  for  layers  of  constant  thickness,  e.g.. 


for  k = 1 and  3,  where  the  limits  k ± 1 denote  constant  z-levels  bear- 
ing those  subscripts,  and  where  h^  is  constant.  The  vertically  averaged 
pressure  is  denoted  by  tt^. 

With  these  definitions  and  approximations,  the  governing  momentum, 
heat,  and  continuity  equations,  vertically  integrated  over  layer  k,  take 
the  following  form: 
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where  T and  are  tctal  vertical  transfer  terms  defined  by 
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lue  ’•ight  sides  of  Eqs.  (3.4)  and  (3.5)  are  evaluated  at  levels  0,  2, 
and  4.  The  advectivt.  fluxes  are  absent  in  the  sea-surface  boundary 
conditions  (2.7)  and  (2.8)  because  there  is  no  advecti ve  flux  (in  a 
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Lagrangian  sense)  across  a fluid  material  surface.  In  the  present 
approximation,  surface  conditions  are  applied  at  mean  sea  level,  and 
continuity  of  fluxes  of  momentum  and  heat  means  that  vertical  advec- 
tions  must  be  included.  In  place  of  Eqs.  (2.7)  and  (2.8)  we  have 

lo  = Is(X*  ° (3.6) 

Qq  = Qs(*.  CP,  t)  (3.7) 

No  modification  is  required  at  the  bottom,  where  Eqs.  (2.9)  through 
k k 

(2.11)  give  and  = 0.  Formulated  in  the  model  is  the 

assumption 


3b  = KH3  0.8) 

where  K is  a constant  bottom-friction  coefficient. 

Evaluation  of  Eqs.  (3.4)  and  (3.5)  at  the  intermediate  level  2 
requires  additional  assumptions.  For  the  advection  terms,  w ^ is  eval- 
uated from  (3.3),  and  v^  and  are  assumed  to  be  simple  two-point 
averages  of  the  torm 

t2  = J <e!  + e3)  (3.9) 

with  a similar  expression  relating  v?  to  and  U^.  (See  Fig.  1.) 

The  above  choice  is  inspired  by  conservation  of  quadratic  quant- 
ities, and  hence  stability  of  the  numerical  scheme,  although  not  in  an 
exact  sense  for  a (divergent)  model  with  a free  upper  surface.  Assuming 
comparable  treatment  of  the  horizontal  finite  differencing,  the  square 
of  a scalar  such  as  T is  conserved  in  the  sense  of  Bryan  (1966)  if  aver- 
ages of  the  form  of  (3.9)  are  used,  provided  also  that  the  normal  com- 
ponent of  velocity  is  zero  at  all  boundaries.  Conservation  of  I2,  U2 , 

2 

and  V is  only  approximate  in  the  present  model  because  the  vertical 
velocity  is  not  zero  at  the  upper  "boundary"  at  mean  sea  level.  Thus, 
for  example, 
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dt  fffT  dxd ydz 


mean  sea 
level 


^here  volume  integration  is  over  the  entire  space-fixed  volume  bounded 
above  by  mean  sea  level.  Conservation  of  T,  U,  and  V,  and  their  squares, 
in  the  present  model  is  in  the  more  limited  sense  that  the  usual  con- 
servation would  hold  if  the  vertical  velocity  were  zero  at  the  upper- 
most level. 

The  vertical  derivatives  in  Eqs.  (3.4)  and  (3.5)  are  evaluated  as 
the  average  of  the  respective  derivatives  occurring  just  above  and  just 
below  level  2.  This  requires  assumptions  on  the  more  detailed  vertical 
distributions  of  v and  T.  Assumed  in  the  model  are  the  piecewise-linear 
distributions  shown  in  Fig.  1.  The  temperature  T^  is  set  equal  to  T^ 
as  a crude  representation  of  the  observed  mixed  upper  layer  in  the 
oceans.  T^  is  set  equal  to  T^  to  be  consistent  with  zero  conduction 
into  the  bottom.  The  temperature  derivative  at  level  2 is 


(3.10) 


in  which  T2  is  given  by  Eq . (3.9).  With  the  assumed  vertical  distribu- 
tion, the  other  temperatures  can  be  related  to  the  averages.  One  finds 


T1  = r (70,  - e,) 


(3.11) 


= T (79,  - 0.) 


(3.12) 


As  stated  previously,  h^  is  assumed  greater  than  the  depth  of  pene- 
tration of  appreciable  solar  radiation.  Therefore,  q = 0 at  level  2. 
This  completes  the  specification  of  Eq.  (3.5)  for  the  vertical 
differencing. 

For  the  velocity  it  does  not  seem  clear  that  a "mixed  layer"  of 
velocity  would  necessarily  be  appropriate.  (Velocities  in  the  general 
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circulation  are  not  readily  subject  to  direct  measurement  except  in 
regions  of  strong  and  persistent  currents.)  Implicit  in  the  model  is 
the  assumption  that  h^  is  much  larger  chan  the  thickness  of  the  upper 
frictional  layer.  Therefore,  Ekman-layer  dynamics  are  not  resolved 
by  the  model.  Perhaps  the  simplest  assumption,  and  the  one  used  here, 
is  to  let  Vq  be  a linear  extrapolation  from  v:  and  v2 . (This  assump- 
tion can  be  easily  modified  in  the  model.)  Again,  v^  is  set  equal  to 
Y 3 m°del  a no-stress  distribution  in  the  lowermost  portion,  and  an 
equation  of  the  form  of  (3.10)  is  assumed  for  the  velocity  derivative 
at  level  2.  The  velocity  is  again  a two-point  average  of  the  form 
of  (3.9),  and  one  finds 


~1  = Hi 


(3.13) 


(3.14) 


This  completes  the  specification  of  (3.4)  for  the  vertical  differencing. 

The  layer  equations  are  completed  with  the  following,  the  first 
two  of  which  are  the  approximations  of  Eqs.  (2.13)  and  (2.14)  when  the 
upoer  limits  of  integration  are  replaced  with  zero. 
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The  above-mentioned  first-order  approximation  is  seen  to  amount  to  a 
linearization.  Equation  (3.15),  together  with  the  continuity  equation 
and  the  condition  w^  = 0 , give  3£/3t  = w^,  which  ir  a linearization  of 
(2.6). 

Equation  (3.16)  is  the  vertical  average  of  (2.14)  over  layer  1,  in 
which  the  upper  limit  of  integration  £ is  replaced  with  zero.  The 
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dynamically  unimportant  term  p^z  integrates  to  a constant,  which  is 
omitted.  The  first  two  terms  on  the  right  of  (3.16)  are  the  contribu- 
tions to  pressure  which  can  occur  in  the  simple  case  of  a homogeneous 
fluid.  Of  these,  the  term  p is  probably  negligible  in  most  simu  itions 

Si 

of  the  general  circulation  (see  Gates  1966)  but  is  retained  here  for 
slightly  greater  generality.  The  last  terms  are  the  baroclinic  contri- 
butions to  pressure.  The  particular  form,  e.g.,  the  rational  factor 
8/18,  resulted  from  the  assumed  piecewise-linear  distribution  of  tem- 
perature. For  reasonable  values  of  temperature,  the  baroclinic  part 
of  (3.16)  is  nearly  equal  to  - y p^gah^,  which  would  have  resulted 
if  temperature  had  been  assumed  vertically  uniform  in  the  top  layer. 

Equation  (3.17)  did  not  result  from  vertical  integration  of  (2.14) 
but,  rather,  is  inspired  by  (approximate)  conservation  of  total  energy. 
Equation  (3.17)  can  be  rewritten  as 

, ' ?■  - • PJ^T2  <3-18> 

2 hl  + 2 h3 

which  is  a finite-difference  approximation  to  the  hydrostatic  equati'n 
evaluated  at  level  2.  It  can  be  shown  (e.g.,  Haney  1971a)  that  Eq . 
(3.18)  represents  one  of  several  (sufficient)  conditions  on  the  ensuing 
numerical  scheme  to  conserve  total  energy  the  vertical  velocity  were 
zero  at  mean  sea  level. 

B.  HORIZONTAL  AND  TIME  DIFFERENCING 

The  foregoing  system  of  equations  is  written  in  finite  difference 
form  on  a horizontally  staggered  grid.  The  grid  is  similar  to  that 
used  by  Bryan  (1969a)  and  others,  except  possibly  for  the  location  of 
coastal  boundaries.  The  basic  unit,  shown  in  Fig.  2,  for  constructing 
the  grid  consists  of  two  temperature  elements  stacked  one  over  the 
other.  Average  temperatures  are  stored  at  the  center  of  each  of  the 
two  elements.  The  average  velocities  are  stored  at  the  corners  of  the 
elements  at  the  same  levels  as  the  temperatures.  Surface  elevation, 
which  is  independent  of  depth,  is  stored  at  a horizontal  position  co- 
incident with  the  temperature  points.  Vertical  velocities  (not  show.) 
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• — U,V  points 


Fig.  2— Schematic  representation  of  two  temperature  box  elements  stocked  one 
over  the  other  This  double  element  is  the  basic  building  block  of  the  two-level 
ocean  model.  Average  temperatures  Qv , shown  in  Fig.  1,  ore  stored  ot  the  center 
eoch  element  on  levels  1 ond  3.  Average  velocities  U|<  ore  stored  ot  the 
corners  of  the  boxes.  Layer  thicknesses  ore  shown  equal^to  each  other  for  clority. 
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are  evaluated  diagnostically  at  levels  2 and  0.  As  described  below, 
this  is  done  separately  for  U,  V points  and  for  0 points. 

The  grid  is  regular  in  a spherical  coordinate  system,  i.e.,  AX 
and  A(p  are  chosen  constant.  The  east-west  distance  between  adjacent 
velocity  points,  and  between  adjacent  temperature  points,  is 
a cos  cpAX,  which  decreases  toward  the  poles.  The  respective  north- 
south  dimensions  are  uniformly  aA(p. 

Horizontal  velocities  are  stored  at  integer-index  values  of  lati- 
tude and  longitude,  in  which  indices  are  co anted  positive  eastward  and 
northward  with  respect  to  reference  values  X^,cp^,  located  at  the  extreme 
southwest  corner: 

>>i  = + (i  - 1)AX,  i = 1,  2,  ...  , M + 1 

<Pj  = ^ + <J  " 1>A<P»  j = 1,  2,  ...  , N + 1 

The  temperatures  and  surface  elevation  are  stored  at  half-integer  index 
locations  given  by 

X .+l  * X.  ± y AX 

i±-J  l 2 

<Tj±.3  = cpj  ± J A(p 

In  general  there  are  (M  + 1)  * (N  + 1)  x 2 velocity  points  and  M x N x 2 
temperature  points,  or,  equivalently,  M * N of  the  double  elements  of 
Fig.  2.  Some  elements  will  denote  artificial  land  points,  to  be  dis- 
cussed presently. 

The  simplest  case  of  a rectangular  basin  is  shown  schematically 
in  Fig.  3.  Indicated  there  by  the  dashed  lines  is  a typical  double 
element  for  temperature  as  viewed  from  above.  The  basin  in  this  case 
is  bounded  by  single  rows  and  columns  of  artificial  temperature  elements. 

The  following  standard  notation  is  adopted  for  divided  differences 
and  two-point  averages: 


O X- 


• = u , v (integer  index  ) points 

0=  (half- integer  index)  points 


:'3-  3 W of  the 

of  Fig.  2.  Showo  IS  the  simple  cose  of  o °"9U  , . denoted  b the  doshed 

"T'l„,il  I ='  111 ' ' si  ’* 
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^ith  similar  expressions  with  respect  to  latitude  and  time.  Thus,  for 
instance,  a central  difference  in  time  evaluated  at  time  t with  incre- 
ment 2 t is  5tfi (t)  . A similar  notation  is  used  for  vertical  divided 
differences  and  averages;  however,  a subscript  k is  used  to  distinguish 
between  th"  upper  and  lower  elements  of  Fig.  2.  For  instance, 


5 w 
2 


k 


w 


k+1 


in  which  k assumes  the  values  1 and  3. 

The  Laplacian  in  finite  difference  form  will  be  the  same  when 
applied  to  both  temperature  and  velocity: 
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a cos  cp 
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cos  cp 
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6 (cos  cr  6 0) 


These  quantities  require  no  averaging  of  the  dependent  variables  at 
intermediate  points,  wherear  the  advections  will  require  averaging. 

The  advection  operator  used  at  temperature  points  differs  from 
that  used  at  velocity  points.  Equations  (3.2),  (3.3),  and  (3.15) 
applied  to  the  elements  of  Fig.  2 are  given  the  following  finite  dif- 
ference forms: 
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(3.19) 
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where 
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and  where  7 


is  obtained  from  (3.22)  by  setting  0 


1 in  the  latter, 


The  form  of  (3.22)  is  quadrati  c-conserving  in  the  sense  of  Bryan 
(1966).  The  horizontal  fluxes  of  G and  of  0 across  vertical  faces  of 
all  temperature  elements  are  equal  and  opposite.  For  instance,  the 
eastward  flux  of  9 across  each  of  the  easternmost  faces  of  the  double 


element  shown  in  Fig.  2 involves  the  product  of  iP,  which  is  the  north- 
south  average  of  the  eastward  velocities  at  the  northeast  and  southeast 
corners,  and  0 , uhich  is  the  east-west  average  of  the  respective  tem- 
peratures of  the  element  shown  and  those  of  the  element  immediately  to 
the  east.  This  flux  is  exactly  equal  (but  of  opposite  sign)  to  that 
entering  the  western  face  of  the  neighboring  element  as  evaluated  by 
(3.22),  For  the  vertical  advective  fluxes  the  quantity  Q*  in  Eq . (3.19) 
contains  the  term  w^  = w^G^  + Q^)/2,  where  w2  is  evaluated  from 
(3.20).  The  vertical  advective  flux  between  the  two  elements  is  again 
equal  but  of  opposite  sign,  and  cancels  out  in  the  summation  process. 

A corresponding  set  of  equations,  but  with  a different  advection 
operator,  is  applied  to  all  interior  velocity  points.  A typical  velocity 
element  viewed  from  above  is  indicated  by  the  dotted  lines  in  Fig.  3. 

Such  elements  are  similar  to  that  shown  in  Fig.  2 except  that  and  U 
are  situated  along  the  center  line  and  C,  0^  and  03  are  at  the  corners. 

Equations  (3.1)  and  (3.3)  are  applied  to  these  elements  in  the  following 
finite  difference  form: 


- V • (UU)  + A A2u  + i-  { T* 
d ~~  k m ~k  p z~k 


r k x u,  - I — i siM 

k~  P^a  ^cos  cp  X k’  <p\  / 


(3.23) 
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6Z”k  ' - Vb  • "k  (3.24) 

where 

7b  * (^k  = a~^s  (p{6X(^V}  + ^[(GQS  cp  V)^P]}k  (3.25) 

and  where  V • Uk  is  obtained  from  Eq . (3.25)  by  replacing  the  vector 
quantities  U and  on  the  right  side  with  one.  The  pressure  dif- 
ferences in  Eq.  (3.23)  are  evaluated  using  (3.16)  and  (3.17). 

The  advection  operator  is  quadratic  conserving  in  and  V . 
Again,  the  vertical  advection  is  w^  - w^  + U3)/2,  which  is  con- 
tained in  the  term  Tj  of  (3.23),  and  is  of  a quadratic  conserving  form. 
Here,  w2  is  calculated  from  (3.24)  instead  of  (3.20). 

At  the  beginning  of  an  integration,  and  for  subsequent  restarts, 
a single  forward  time  difference  is  used  in  place  of  the  central  time 
differences  in  Eqs . (3.19),  (3.21),  and  (3.23).  Central  time  differ- 
ences are  used  for  the  remainder  except  that  the  diffusion  terms  are 
lagged  one  time  step  for  stability. 
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IV.  GENERAL  GRID  SETUP  AND  BOUNDARY  CONDITIONS 

A,  TREATMENT  OF  IRREGULAR  BOUNDARIES 

The  basic  idea  is  to  fill  all  water  regions  of  an  ocean  basin  with 
the  temperature  double-elements  of  Fig.  2 in  such  a way  that  velocity 
points  lie  on  the  boundaries.  The  (horizontal)  velocities  are  held 
identically  zero  there.  Land  regions  are  filled  with  similar  but  arti- 
ficial elements  which,  except  for  land-coastal  points,  are  never  used 
in  the  calculations.  At  all  land-coastal  points,  i.e. , t 'mperature- 
box  elements  immediately  adjacent  to  water  regions,  the  values  of 
temperature  are  adjusted  to  satisfy  the  insulating  condition  at  coasts 
(meaning  zero  conduction  there  because  norizor.tal  advections  are  zero 
from  the  velocity  condition).  The  approach  appe  rs  similar  to  that 
used  in  the  Indian  Ocean  model  of  Cox  (1970). 

The  insulating  condition  is  satisfied  approximately  in  the  present 
model  by  setting  land-coastal  values  of  temperature  equal  to  the  average 
of  values  for  adjacent  sea-coastal  points.  This  is  done  initially  and 
at  the  end  of  every  subsequent  time  step.  The  method  is  indicated  in 
Fig.  4 for  the  case  of  an  island.  The  temperatures  0^  and  0^  at  the 
land  point  labeled  8 are  set  equal  to  the  corresponding  temperatures 
at  the  sea  point  immediately  to  the  left,  which,  in  this  case,  exactly 
satisfies  the  insulating  condition.  However,  the  temperatures  at  the 
land  point  labeled  q + 1 are  set  equal  to  the  averages  ^ and  03  re- 
spectively) of  the  temperatures  at  sea-points  labeled  9 and  1,  which 
is  only  an  approximation  to  the  insulating  condition.  An  extreme  case 
is  that  of  an  island  consisting  of  a single  element  for  which  the 
temperatures  and  0^  are  the  respective  averages  of  those  for  the 
four  adjacent  sea  points.  Outer  lateral  boundaries  of  the  bounded 
basin,  not  shown  in  Fig.  2,  are  treated  in  a similar  fashion. 

It  is  also  possible  to  formulate  the  model  so  as  to  satisfy  exactly 
the  insulating  condition  at  coast  lines,  but  the  method  is  more  complex 
than  seems  justified  in  a model  that  neglects  continental  shelf  topo- 
graphy. In  the  exact  method  a significant  number  of  land-coastal  tem- 
peratures would  have  to  be  changed  twice  (some  as  many  as  four  times) 
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Fig.  4 — Schematic  indicoting  the  method  used  in  constructing  the  more  generol 
grid  ond  in  sotisfying  loterol  boundory  conditions  for  irregulor  geometries.  The 
ocean  is  filled  with  the  temperoture  elements  of  Fig,  2.  Adjacent  lond  oreas 
(e.g.,  the  islond  shown)  ore  filled  with  similor  but  ortificial  elements.  Vel- 
ocities U,  V ore  maintained  identically  zero  within  oil  lond  oreas  ond  on  the 
boundaries.  Insuloting  conditions  ot  coastlines  are  satisfied  approximately  by 
setting  lond-coastol  volues  of  temperoture  equal  to  the  overage  of  volues  for 
adjacent  sea-coastal  points.  See  text  for  explanotion  of  numbers. 
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during  the  calculations  for  each  time  step.  In  any  event,  the  bound- 
aries for  real  ocean  circulations  are  not  well  defined  her<  and  are  not 
precisely  insulated.  (Tn  the  case  of  the  Pacific  model  reported  on 
below,  the  "boundaries"  were  placed  at  the  100-m  depth  contour.)  On 
balance,  the  approximate  method  appears  better  because  of  its  simplicity. 

B.  THE  GEOGRAPHY  ARRAY 

The  general  treatment  of  boundary  conditions  is  perhaps  best  de- 
scribed in  terras  of  the  geography  array  which  is  input  to  the  model. 

This  is  an  M * N array  corresponding  to  the  total  number  of  temperature 
double-elements  of  Fig.  2 (cf.  the  simple  M * N array  of  Fig.  3).  Each 
element  of  the  input  array  contains  one  of  three  simple  pieces  of  infor- 
mation: zero  (blank)  if  it  is  a regular  sea  point,  "I"  i . it  is  a sea 

point  covered  by  sea  ice,  or  "L"  if  it  is  a land  point.  The  program 
then  processes  this  information  and  sets  up  a slightly  mere  elaborate 
M x N array  which  controls  the  treatment  of  each  point  at  every  time 
step.  The  initial  setup  goes  somewhat  as  follows. 

For  all  land  points  the  associated  velocities,  eight  for  each 
element  shown  in  Fig.  2,  are  set  to  zero  and  tagged  so  that  they  are 
never  changed  during  the  course  of  the  integration.  Land-temperature 
points  away  from  the  coast  are  similarly  tagged  and  never  used  in  the 
calculations.  All  land-coastal  points  are  given  a special  number, 
stored  in  the  geography  array,  according  to  the  type  of  point.  This 
number  tells  the  program  how  to  do  the  temperature-averaging,  initially 
and  at  the  end  of  each  time  step.  For  instance,  all  land-coastal  points 
such  as  the  one  labeled  2 + 3 in  Fig.  4 are  given  the  same  identifying 
number.  After  every  time  step  the  temperatures  at  all  such  2+3  points 
are  set  equal  to  the  averages  of  the  respective  temperatures  at  the  sea- 
points  immediately  to  the  east  and  to  the  north. 

All  sea-ice  points  retain  a single  identifying  number  and  are 
treated  the  same  way.  The  temperatures  at  both  levels  are  set  and  held 
fixed  at  -1.9  C,  roughly  the  (laboratory)  equilibrium  temperature  of 
ice  and  saline  water  of  35  /QO  salinity.  The  temperature  prediction 
equation  is  ignored  at  these  points,  as  is  the  (input)  surface  heat 
flux  Qg.  However,  the  surface  wind  stress  Tg  together  with  Eqs.  (3.23) 
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and  (3.24)  are  used  to  predict  the  velocities  at  these  points.  Sea 
ice  is  an  exceedingly  inhomogeneous  and  anisotropic  substance  contain- 
ing numerous  cracks  and  leads,  and  is  in  nearly  constant  motion  from 
surface  winds.  The  model  represents  sea  ice  as  a massless,  completely 
extensible  and  compressible  membrane  through  which  wind  stress  is  trans- 
mitted without  attenuation  to  the  nearly  homogeneous  water  below.  This 
is  probably  one  of  the  cruder  features  of  the  model. 

finally,  if  the  geography  array  indicates  a regular  sea  point,  the 
complete  set  of  prediction  equations  is  employed,  including  (input) 
surface  heating  and  wind  stress.  All  sea  points  adjacent  to  land, 
both  velocity  and  temperature  points,  are  treated  just  like  the  corre- 
sponding points  in  the  open  ocean. 
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V.  INITIALIZATION 


A.  THE  PROBLEM  AND  SUMMARY  OF  THE  METHOD 

The  purpose  of  the  model  is  similar  to  that  of  atmospheric  pre- 
diction models.  For  given  initial  distributions  of  the  dependent 
variables  (in  this  case  temperature,  surface  height,  and  velocity)  the 
model  is  to  make  an  oceanic  "forecast"  of  the  dependent  variables, 
subject  in  this  case  to  the  prescribed  time-evolution  of  surface  heat- 
ing and  wind  stress.  However,  as  is  well  known  in  atmospheric  fore- 
casting, the  forecast  problem,  particularly  with  the  unfiltered  primitive 
equations,  introduces  a secondary  problem  of  initialization  which  in 
some  respects  is  more  difficult  to  solve.  Initial  temperature  distri- 
butions can  be  prescribed  reasonably  well  from  observation.  However, 
except  for  isolated  regions  of  strong  and  persistent  currents,  the 
velocity  field  In  the  general  circulation  of  the  oceans  is  not  a 
directly  observable  quantity.  It  is  perhaps  less  obvious  but  nonethe- 
less true  that  the  small  spatial  variations  in  pressure  (or  surface 
elevation)  that  are  important  in  the  oceanic  circulation  are  also  dif- 
ficult to  measure.  Related  to  the  fundamental  observational  difficulty 
is  the  somewhat  more  practical  problem  (from  a computational  point  of 
view)  of  smoothing  or  filtering  the  input  data.  The  governing  equa- 
tions of  the  present  model  do  not  filter  out  surface  gravity  waves. 
Arbitrarily  specified  initial  conditions  (even  if  known)  can  generate 
unrealistically  large  surface  waves  whose  amplitudes  may  be  comparable 
to  the  ocean  depth,  for  instance.  These  constitute  oceanographic  noise 
which,  to  avoid  contaminating  th  : ensuing  prediction,  must  be  removed 
by  an  initialization  scheme. 

The  problem  of  initialization  for  atmospheric  prediction  was 
apparently  first  investigated  by  Charney  (1955).  Many  schemes  have 
evolved  Si.nce  that  time;  for  example,  the  global  atmospheric  initiali- 
zation scheme  of  Houghton  and  Washington  (1969).  However,  no  comparable 
scheme  for  ocean  models  appears  to  have  been  published.  Initialization 
was  not  considered  in  the  published  work  of  Crowley  (1968). 
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The  scheme  worked  out  here  is  an  extension  of  the  diagnostic 
method  of  calculating  currents  from  observed  density  fields  as  form- 
ulated by  Sarkisyan  and  coworkers  (e.g.,  Sarkisyan  and  Ivanov  1971). 
First,  an  initial  state  for  surface  elevation  and  velocity  is  deter- 
mined for  a given  fixed  temperature  distribution.  Here,  use  is  made 
of  two  fundamental  characteristics  of  the  circulation  which  distinguish 
it  from  gravity-wave  behavior,  as  pointed  out  by  Charney  (e.g.,  Charney 
1955)  : its  quasi-r.ondivergence  (in  the  horizontal)  over  large  horizontal 

distances,  and  its  quasi-geost rophy  on  the  rotating  earth.  The  (geo- 
strophic)  thermal  wind  (or  thermal  current)  equations  appear  to  be 
sufficient  to  introduce  a baroclinic  component  into  the  flow  field 
without  generating  oceanographic  noise  (verified  empirically). 

The  second  part,  which  is  analogous  to  the  diagnostic  method  of 
Sarkisyan,  consists  of  a time  integration  employing  the  more  complete 
set  of  momentum  equations.  This  allows  for  "barostropic  spinup,  or 
geostrophic  adjustment  of  the  pressure  and  velocity  fields  to  the  pre- 
scribed density  field  and  boundary  conditions.  Temperature  and  wind 
stress  are  held  fixed  and  heating  is  ignored  throughout  the  initializa- 
tion process . 

Thirdly,  the  fields  are  (or  may  be)  time-averaged  to  produce  an 
initialization  stats  which  is  quasi-steady  for  the  given  fixed  dis- 
tributions of  temperature  and  wind  stress.  By  this  is  meant  that  the 
time-averaged  flow  is  in  dynamical  equilibrium,  or  nearly  so.  In  gen- 
eral, there  will  not  be  thermal  equilibrium,  as  indeed  is  probably  the 
case  in  the  real  ocean.  When,  after  initialization,  heating  and  the 
temperature-prediction  equation  are  included,  the  temperature  will  gen- 
erally change.  This,  in  turn,  causes  the  velocity  to  change. 

B.  INITIAL  STATE 

The  initial  flow  is  to  be  nondivergent  and  geostrophic,  i.e., 

2£.=  -y»Hv=0 

at 

fk  x v = - ~~  Vp 
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whore  the  overbars  denote  vertical  averages  over  the  total  depth  H. 

It  is  sufticient  here  to  take  v = 0 and  p = const  as  initial  conditions. 
The  latter,  in  terms  of  the  layer-average  pressures,  is 

h^i  + h ^ = const  (5.1) 

which,  when  Eqs.  (3.16)  and  (3.17)  are  substituted,  gives  a relation- 
ship between  the  initial  surface  elevation  and  the  initial  prescribed 
temperature  distribution.  The  above  constant  is  adjusted  to  give  zero 
mean  surface  elevation. 

The  condition  v = 0 gives  one  vector  relationship  between  the 
average  velocities  for  each  of  the  two  layers: 

hlUl  + h3U3  = ° (5‘2) 

A second  relationship  comes  from  the  thermal  wind  (or  thermal  current) 
equations.  These  are  the  vertically  differentiated  geostrophic  equa- 
tions in  which  pressure  is  replaced  with  temperature  by  means  of  the 
hydrostatic  equation  and  the  equation  of  state.  In  the  present  model 
the  vector  thermal  rfind  equation  takes  the  finite-difference  form 

fk  x (u  - U3)  * - j gHaT2  (5.3) 

where  T2  is  the  arithmetic  average  of  9^  and  9^. 

Equations  (5.1)  through  (5.3)  give  the  initial  fields  of  C,  Up 
and  for  a specified  initial  temperature  field.  However,  such  fields 
must  be  modified  slightly  to  satisfy  the  boundary  conditions.  As  de- 
scribed in  the  previous  section,  land-coastal  values  of  temperature  are 
adjusted  to  satisfy  (approximately)  the  insulating  condition  at  coast- 
lines, and  coastal  values  of  and  are  set  to  zero. 

C.  BAROTROPIC  SPINUP  AND  TIME  AVERAGING 

The  above  initial  state  is  found  a posterior ” to  be  sufficient  to 
introduce  a baroclinic  component  in  the  flow  field  without  generating 
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un realistic  gravity-inertial  oscillations.  To  complete  the  initiali- 
zation we  must  adequately  simulate  the  geostrophic  adjustment  which 
can  take  place  in  a barotropic  (or  homogeneous)  fluid  giving  a component 
of  velocity  independent  of  depth  (and  v 4 0 in  general). 

A description  of  the  spinup  of  a homogeneous  ocean  in  a shallow 
basin  of  planetary  scale  is  given  by  Gates  (1968).  We  merely  point  out 
here  that  characteristic  features  of  the  adjustment  process  are  the 
westward  intensification  of  currents  (the  "Gulf  Stream")  accompanied  by 
westward-propagating  planetary  or  Rossby  waves.  Spinup  times  (measured 
in  terms  of  the  first  peak  in  kinetic  energy)  are  of  the  order  of  a few 
tens  of  days  to  a few  months.  Surface  gravity  waves  of  more  moderate 
amplitude  are  in  general  present.  Transit  times  for  the  principal 
gravity  wave  modes  are  of  the  order  of  2L/(gH)1/2  ~a  few  days  to  a 
few  weeks,  where  L is  a typical  horizontal  dimension  of  the  basin. 

The  purpose  here  is  not  to  study  transient  phenomena,  but  rather 
to  generate  initial  fields  of  velocity  and  surface  elevation  which  are 
quasi-steady  with  respect  to  given  initial  fields  of  temperature  and 
wind  stress.  This  is  achieved  in  practice  by  employing  the  foregoing 
initial  state  and  performing  a time-integration  of  the  primitive  equa- 
tion model  for  a period  spanning  about  twice  the  spinup  time,  while 
holding  the  temperature  and  wind  stress  fixed.  At  this  point  the 
energies  are  inspected  visually;  they  typically  show  peaks  and  troughs 
which,  beyond  spinup,  represent  oscillations  between  kinetic  and  po- 
tential energies  as  the  sura  K + P approaches  a constant.  For  simple 
rectangular  basins,  a dominant  period  — 2L/(gH)^^2  usually  shows  up 
clearly  in  such  oscillations , provided  bottom  friction  is  small.  For 
more  complex  basin  geometries  there  may  be  no  single  well-defined 
period.  There  mav  be  a composite  of  many  periods,  or  there  may  be 
none  discernible  (the  latter,  e.g.,  for  a highly  viscous  and  damped 
case).  If  peaks  and  troughs  are  discernible  in  the  energies,  the  model 
is  restarted  and  time-averaging  performed  between  an  integral  number  of 
chosen  peaks  (say)  to  remove  the  oscillations.  The  time -averaged  fields 
are  then  substituted  back  and  the  model  restarted.  If  the  energies 
continue  to  show  some  appreciable  oscillations,  the  process  is  repeated 
until  a quasi-steady  state  is  achieved  for  the  prescribed  temperature 
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and  wind  stress  distributions.  This  represents  the  initialization 
state , from  which  predictions  of  temperature  and  velocity  are  subse- 
quently made  with  the  more  complete  set  of  equations. 

The  above  procedure  has  been  applied  to  a number  of  cases  ranging 
from  a simple  rectangular  basin  to  the  more  complex  Pacific  model  re- 
ported on  ■ elow.  Typical  results  in  terms  of  the  energies  for  cases 
i.  solving  simpler  basin  geometries,  and  without  appreciable  frictional 
damping  of  oscillations,  are  shown  schematically  in  Fig.  5.  The  geo- 
strophic  and  nondiveigent  initial  state  is  input  at  time  zero.  Reason- 
ably well-defined  oscillations  occur  between  kinetic  and  potential 
energies  after  spinup.  For  a simple  shallow  basin  with  horizontal 
scale  L ~ 5000  km,  a time  span  for  averaging  (B  - C)  of  10  to  15  days 
is  usually  sufficient.  Typically,  only  one  time-averaging  period  was 
required  to  achieve  sufficiently  steady  flow.  The  initialization  state 
is  at  time  C and  beyond,  so  long  as  temperature  and  wind  stress  remain 
fixed. 

Computer  time  requirements  for  this  model  in  the  complete  predic- 

-4 

tive  phase  on  an  IBM  360/91  computer  are  less  than  2 x 10  sec  per 
time  step  per  double  temperature  element  of  Fig.  2.  This  amounts  to 
approximately  4 minutes  of  machine  time  per  model-day  for  the  Pacific 
model  reported  on  in  Sec.  VI.  This  value  compares  favorably  with  the 
approximate  value  of  14  minutes  per  model-day  required  on  the  same 
machine  by  an  improved,  faster  version  of  the  Mintz-Arakawa  global  atmos- 
pheric model.  Further  discussion  of  this  point  is  given  in  Sec.  VII. 
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VI.  THE  PACIFIC  MODEL 


! 


The  basic  model  for  an  arbitrary  ocean  basin  was  described  in  the 
foregoing  sections.  The  purpose  of  this  section  is  to  present  results 
from  a demonstration  run  for  the  case  of  the  Pacific. 


A.  INPUT  DATA 

The  basin  extends  from  80°S  to  64°N.  Artificial  meridional  bar- 
riers are  placed  at  105  E (just  west  of  Australia)  and  across  Drake's 
Passage  at  67  W.  The  geography  array  described  in  Sec.  IV  was  set  up 
based  on  the  data  of  Smith  et  al.  (1966),  who  compiled  global  distribu- 
tions of  land  elevation  and  sea  depth  on  approximately  a one-degree  grid. 
Because  the  present  model  neglects  all  topographic  effects,  I chose 
(somewhat  arbitrarily)  the  100-m  depth  contour  as  the  land-sea  boundary. 
Grid  increments  are  two  degrees  of  latitude  ar.d  longitude.  There  re- 
sulted 73  y 96  of  the  double  elements  of  Fig.  2,  including  all  artifi- 
cial land  elements.  The  array  includes  part  of  the  western  North 
Atlantic  which  should  be  ignored  in  the  results.  The  array  was  also 
defined  in  such  a way  that  temperature  points  lie  on  the  equator. 
Therefore,  velocity  points  are  displaced  a distance  A<f>/2  away  from  the 
equator,  which  avoids  singularities  in  the  geostrophic  initial  velocities. 

Fixed  parameters  include  the  following: 


Horizontal  diffusion: 
Vertical  diffusion: 
Bottom  friction: 

Layer  thicknesses: 
Grid  size: 

Time  increment: 


Am  = 5 x 108cmVsec,  A^  = 5 * lC7cm2/sec 

2 2 
v = 10  cm  /sec,  < = 0.3  cm  /sec 

K = 2 * 10  dyne-sec/cm8 

(coefficient  of  2?  10  ^/sec) 

h1  = 100  m,  h3  = 200  m 

AA  = A$>  = 2° 

At  = 7.5  min 


The  above  value  of  A^  is  controlled  somewhat  by  the  grid  size,  and 
places  the  solutions  in  the  viscous  regime  studied  by  Munk  (1950).  The 
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total  depth 
level  model 
(Cromwell) 


of  300  m is  of  the  order  of  magnitude  required  for  a two 

to  resolve  the  strong  but  relatively  shallow  equatorial 
undercurrent. 


1 he  above  values  of  v and  < were  chosen  also  in  order  to  model 
equatorial  regions  as  accurately  as  possible.  Arthur  (1960)  and  Charney 
(I960)  have  recommended  values  of  v of  order  10  cm2 /sec,  while  Wyrtki 
and  Bennett  (1963)  and  Montgomery  (as  reported  by  Wyrtki  and  Bennett) 
have  recommended  values  of  < of  a few  tenths.  Such  relatively  small 
values  of  k-  are  also  consistent  with  thermocline  theory  for  middle- 
latitude  regions  of  the  open  ocean  (e.g.,  Alexander  1971).  (Wyrtki  and 
Bennett  recommend  v - 3 cm2/sec:  however,  in  the  present  model  v is 
more  like  the  vertically  averaged  quantity  deduced  by  Arthur.) 

The  inclusion  of  bottom  friction  was  motivated  by  a desire  to  re- 
duce the  number  of  calculations  required  for  initialization  in  a demon- 
stration run.  Gates  (1972a)  has  shown  that  bottom  friction  significantly 
dampens  the  transients  generated  in  the  spinup  of  a homogeneous  ocean. 
This  led  me  to  include  bottom  friction  under  the  presumptions,  first, 
that  the  (reduced)  energy  peaks  will  be  reached  sooner,  and,  secondly, 
that  time-averaging  will  be  avoided  because  the  transients  will  have 
almost  disappeared  by  the  end  of  spinup.  By  a simple  scale  analysis, 
the  above  value  of  K gives  bottom  friction  forces  of  the  same  order  of 
magnitude  as  horizontal  friction  forces. 

Other  input  data  include  initial  temperature  fields  and  sea-surface 
distributions  of  wind  stress  and  heating.  These  are  based  on  January 
(or  northern  winter)  conditions,  motivated  by  comparable  control  and 
experimental  runs  made  with  a two-level  atmospheric  model  (e.g.,  Gates 
1972b,  Warshaw  and  Rapp  1973).  Initial  temperature  fields  are  shown 
in  Fig.  6 based  on  January  normal  sea-surface  temperature  data  obtained 
from  U.S.  Fleet  Numerical  Weather  Central  (FNWC)  for  the  northern 
hemisphere,  and  National  Center  for  Atmospheric  Research  (NCAR)  for  the 
southern  hemisphere  (as  interpolated  from  updated  global  monthly  normals 
in  preparation  by  Alexander  and  Mobley  1973).  Extrapolation  to  the 
lower  levels  of  the  model  is  based  on  Eq . (3.11),  which  relates  01  and 
3 to  Tj,  and  hence  to  Tq  (cf.  Fig.  1).  a second  relation  between  0 
and  ^ is  based  on  thermocline  theory  in  which  I assumed 


(o ) Vertically  overaged  temperature  for  the  upper  loyer,  with  contour  intervals 
of  2°C  and  the  23-degree  isotherm  dashtd.  The  closed  isotherms  south  of 
Mexico  Yucoton  ore  25°C,  ond  those  northeast  of  New  Guineo  are  27°C 
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63  = Oj  exp(-Cz  /D|sin<()|) 

where  D = 300  m is  a scale-depth,  and  Z3  = hj.  + £ h3  = 200  m is  the 
depth  of  the  lower  level.  With  the  above  scale-depth  it  has  been  shown 
(Alexander  1971)  that  a choice  of  the  constant  C -0.2  can  give  results 
in  reasonable  agreement  with  observation.  Pecause  of  the  singularity 
at^the  equator,  sin*  was  replaced  with  sin  10°  for  latitudes  between 
10  N and  10°S.  This  is  the  reason  for  the  slight  discontinuity  in  the 
ncrth-south  gradient  of  in  Fig.  6b,  noticed  especially  at  10°S . 

Temperature  distributions  appear  reasonable,  except  possibly  for 
the  eastern  equatorial  region  in  the  upper  layer  and  low  latitudes  in 
general  for  the  lower  layer.  (Compare,  e.g.,  Sverdrup  et  al. , 1942, 

Charts  II  and  IV.)  At  the  upper  level,  Fig.  6a,  a tongue  of  relatively 
cold  water  associated  with  the  dashed  isotherm  (~  24°C  at  the  surface 
or  ~23°C  average  for  the  top  layer)  may  not  be  sufficiently  well  de- 
fined along  the  equator  toward  the  east.  One  reason  for  this  may  be 
the  crudeness  of  the  present  grid,  and  the  crudeness  of  the  grids  from 
which  the  normal  temperature  data  were  obtained.  Also,  equatorial 

regions  are  near  the  outer  boundary  of  the  FNVC  northern-hemisphere 
grid. 

Boundary  data  for  the  model  are  shown  in  Fig.  7.  Figure  7a  gives 
th°  wlnd  stress  distribution  according  to  Hellerman  (1967).  Clearly 
evident  are  the  familiar  northeasterly  (southwestward  directed)  trade 
winds  in  low  northern  latitudes.  Winds  along  the  equator  are  generally 
westward  over  the  eastern  two-thirds  or  the  basin.  A maximum  wind 
stress  of  2.G  dynes/cm  occurs  southwest  of  Australia.  We  extrapolated 
Hellerman 's  data  somewhat  into  high  southern  latitudes.  However,  the 
southernmost  two  rows  of  grid  points  (75°S  and  77°S)  contain  zero  values 
for  lack  of  other  information. 

Shown  in  Fig.  7b  is  the  January  normal  distribution  of  surface 
heating  after  Budyko  (1963),  as  interpolated  from  the  4x5  deg  grid 
data  of  Schutz  and  Gates  (1971).  Consistent  with  northern  winter  con- 
ditions, the  North  Pacific  is  generally  being  cooled,  and  the  South 
Pacific  (except  high  southern  latitudes)  generally  heated.  Maximum 
cooling  of  about  780  ly/day  occurs  just  east  of  Japan.  Again,  because 

' 

' 
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data  were  not  available  in  high  southern  latitudes,  the  heating  there 
was  arbitrarily  set  to  zero.  This  may  not  be  important  in  the  model 
because  much  of  the  region  is  covered  by  the  main  ice  pack.  Here,  as 
described  in  Sec.  IV,  the  temperature  is  held  fixed  at  -1.9°C  and  the 
heat  equation  ignored.  (Compare  the  ice-covered  region  in  Fig.  6.) 

B.  RESULTS 

ihe  initial  state  fo?  surface  elevation  and  the  velocities  was 
determined  employing  the  procedure  of  Sec.  V for  the  given  fixed  tem- 
perature distributions  of  Fig.  6.  The  resulting  surface  elevation  is 
shown  in  Fig.  8a.  Note  the  rough  similarity  between  this  pattern  of 
isolines  and  the  pattern  of  isotherms  shown  in  Fig.  6a.  This  results 
from  the  initial  condition  (5.1)  relating  surface  height  to  the  tem- 
perature distribution.  The  initial  velocities,  not  shown,  give  equal 
and  opposite  transports  <n  each  layer  from  the  initial  condition  (5.2). 

Spinup  proceeded  subject  to  the  wind  stress  of  Fig.  7a,  with 
temperature  held  fixed  and  surface  heating  and  the  heat  equation  ignored. 
As  expected,  bottom  friction  was  apparently  sufficient  to  damp  the 
oscillations  generated  during  spinup.  Kinetic  energy  reached  a peak 
near  the  end  of  day  79  and  then  began  to  decline  slightly . The  be- 
havior of  the  energies  near  spinup  is  shown  on  an  expanded  scale  in 
Fig.  9.  The  change  in  total  energy,  kinetic  plus  potential,  is  less 
than  0.1  between  day  74  and  day  90.  Inspection  of  the  fields  at  5-day 
intervals  confirmed  that  the  flow  was  sufficiently  steady  after  about 
day  74  to  represent  an  initialization  state  for  the  prescribed  tempera- 
ture distribution.  (Such  inspection,  taken  by  itself,  is  hardly  proof 
of  steadiness  because  5-day  periodicity,  for  example,  would  not  show 
up.  However,  no  5-day  or  shorter-period  oscillations  are  apparent  in 
Fig.  9.) 

Surface  elevation  at  the  end  of  90  days  is  shown  in  Fig.  8b.  Well- 
definej  subtropical  high-pressure  cells  appear  in  northern  and  southern 
middle-latitudes,  and  the  northern  subpolar  low-pressure  cell  is  prob- 
ably realistic.  Not  much  significance  should  be  attached  to  the  south- 
ern, very  deep,  low-pressure  cell  because  of  the  artificial  meridional 
barriers.  One  encouraging  factor  is  the  more  diffuse  nature  of  the 


denotes  zero.  A maximum 


Boundary  data  for  the  model 


(a)  Initial  surface  elevation  based  on  non-di vergent  and  geostrophic  ‘nitial 
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(b)  Surface  elevation  at  the  end  of  «0  days  (initialization  state)  with 
temperature  held  fixed.  Note  the  smaller  gradients  eost  of  Australia 
^rom£ared^to^those_KW^e  Current. 
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Detailed  behavior  of  energies  near  spinup.  Note  the  expanded  scales. 
Potential  energy  at  day  90  is  0.8%  less  than  its  value  at  day  74. 
Temperature  and  wind  stress  are  held  fixed  throughout  initialization . 
The  more  appreciable  oscillations  depicted  schematically  in  Fig.  5 do 
not  occur  in  this  run  because  of  bottom-friction  damping. 
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high-pressure  cell  east  of  Australia  compared  with  its  counterpart  to 
the  north,  indicating  that  the  Hast  Australia  Current  will  be  weaker 
than  the  Kuroshio,  in  rough  agreement  with  nature.  However,  the  strong 
gradient  at  the  meridional  barrier  west  of  Australia  indicates  a totally 
unrealistic  current  there,  as  expected.  Along  the  equator  the  surface 
elevation  increases  toward  the  west,  consistent  with  the  piling-up  of 
water  under  a westward-directed  wind  stress. 

The  small-scale  diamond  patterns  and  sawtooth  patterns  in  Fig.  8b, 
e.g.,  west  and  east  of  the  Philippines,  respectively,  are  computational 
and  not  physical.  They  represent  slight  solution  separation  that  is 
inherent  in  the  central  space-difference  scheme  used,  apparently  gen- 
erated by  the  computational  boundary  conditions.  Such  separation 
appears  by  the  end  of  day  20  in  tl  * present  integration,  but  does  not 
worsen  significantly  throughout  the  subsequent  calculations. 

The  resulting  initialization  state  for  the  velocities  is  shown  in 
Fig.  10.  Probably  the  most  realistic  region  is  the  North  Pacific.  Tne 
relatively  strong  Kuroshio  current  appears  near  the  Philippines  and 
Japan.  Its  maximum  values  are  nearly  60  cm/sec  at  both  levels;  these 
are  somewhat  on  the  low  side  compared  with  observed  values  of  1 to 
2 m/sec,  but  are  consistent  with  the  crudeness  of  a two-degree  grid 
resolution.  The  Kuroshio  extension,  roughly  the  region  extending  east- 
ward to  the  10  cm/sec  isoiine,  is  reasonably  well  defined,  as  is  its 
further  extension  into  the  North  Pacific  Current,  generally  the  region 
bounded  by  the  dashed  (5  cm/sec)  isoline.  The  Alaska  current  and  the 
southward-flowing  Oyashio  to  the  west  appear  reasonable.  (Compare, 
e.g.,  Sverdrup  et  al.,  1942,  Chart  VII.)  The  southward  drift  in  eastern 
middle  latitudes  might  be  termed  the  California  Current  although  it 
does  not  extend  southward  of  California  in  the  model.  This  may  be  con- 
sistent with  nature  in  the  northern  winter  (e.g.,  the  discussion  in 
Sverdrup  et  al.,  1942,  p.  725,  of  the  northward-flowing  Davidson  Current). 
The  westward-directed  North  Equatorial  Current  is  well  defined  in  lati- 
tudes of  about  15°-20°N  in  the  model. 

The  flow  in  low  latitudes  is  unrealistic  in  several  ways  (partly 
corrected  in  the  ensuing  prediction).  Principally,  the  flow  may  not  be 
sufficiently  east-west,  but  has  significant  poleward  components  in  the 


<"al  ly  - averaged  ve 


(b)  Vertically  overoged  velocity  for  *he  lower  layer.  Maximum  speed  occur* 
•n  the  Kuroshio  and  is  nearly  60  cm/ sec  at  both  levels. 
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Fig.  10  — Initialization  state  for  the  velocity  fields.  The  direction  of  flow 
is  given  by  the  half-arrows,  with  speed  contoured  in  increments 
of  5 cm/sec;  the  5 cm /sec  isoline  is  dashed. 
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upper  layer,  equatorward  in  the  lower  layer,  as  can  be  seen  in  Figs.  10a 
and  b,  respectively.  Slightly  realistic  aspects  are  the  larger  speeds 
depicted  by  the  isolines  near  the  equator,  and  the  tendencies  for  equa- 
torial velocities  to  have  some  westward  component  in  the  upper  layer, 
and  some  eastward  component  in  the  lower  layer.  As  demonstrated  pres- 
ently and  discussed  in  Sec.  VII,  the  unrealistic  meridional  components 
near  the  equator  are  caused  by  an  inaccurate  specification  of  initial 
temperatures  in  low  latitudes. 

The  flow  in  eastern  low  latitudes  is  generally  confusing.  In  low 
southern  latitudes  a northward-flowing  "Peru  Current"  has  left  the 
coast  at  central  Chile.  If  anything,  there  is  southward  flow  (in  the 
upper  layer)  off  Peru,  suggesting  severe  El  Nino  conditions  to  an  extent 
probably  never  observed  (cf.  Sverdrup  et  al.  , 1942,  p.  704). 

Encouraging  features  in  the  South  Pacific  include  the  weaker  and 
more  diffuse  nature  of  the  East  Australia  Current  compared  with  the 
Kuroshio.  However,  the  strong  currents  at  the  meridional  barrier  west 
of  Australia  are  unrealistic  as  expected,  as  is  the  flow  in  high 
southern  latitudes. 

A fifteen-day  prediction  is  now  made  in  the  sense  that  temperature 
is  allowed  to  change  subject  to  the  heat  equation,  and  subject  to  the 
surface  heating  of  Fig.  7b  and  the  (approximately)  insulating  boundary 
conditions.  However,  the  temperature  at  ice-covered  points  continues 
to  be  fixed  at  -1.9°C  as  described  in  Sec.  IV.  The  wind  stress  of 
Fig.  7a  is  also  held  fixed  for  this  experiment.  The  resulting  changes 
that  occur  in  velocity  and  temperature  are  most  pronounced  in  equatorial 
and  eastern  tropical  regions. 

The  temperature  distributions  at  the  end  of  day  105  are  shown  in 
Fig.  11.  The  most  striking  feature  is  the  tongue  of  relatively  cold 
23-deg  water  that  has  developed  in  the  upper  layer  in  the  eastern 
equatorial  region  (Fig.  11a).  The  tongue,  which  is  centered  almost 
exactly  on  the  equator,  appears  to  be  in  rough  agreement  with  Chart  II 
(February  distribution)  of  Sverdrup  et  al.  (1942).  (In  Fig.  11a,  23°C 
corresponds  roughly  to  24-deg  surface  water.)  Temperature  reductions 
in  excess  of  1 C have  occurred  in  the  eastern  equatorial  region.  The 
region  of  27-deg  water  to  the  west  has  been  reduced  slightly  in  size 
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compared  with  the  initialization  state,  and  the  maximum  temperature  is 
now  27. 9°C,  compared  with  28.2°C  before  the  prediction.  This  feature 
may  be  somewhat  less  in  agreement  with  observati>  . 

In  the  lower  layer  (Fig.  lib)  the  situation  is  similar  in  that 
most  of  the  changes  occur  in  low  latitudes.  The  tongue  of  13-deg  water 
toward  the  east  has  increased  in  size  and  has  become  more  nearly  centf  r°d 
on  the  equator.  Toward  the  west  the  region  of  water  with  temperatures 
in  excess  of  15°C  has  become  larger.  (The  region  bounded  by  the  dashed 
line  with  temperatures  less  than  15°C  is  now  smaller.)  However,  a 
cooling  trend  has  occurred  in  the  immediate  vicinity  of  the  equator  as 
evidenced  by  the  westward  extension  of  the  15-deg  water. 

Generally , at  both  levels  the  equatorial  region  has  been  cooled, 
and  at  the  same  time,  the  east-west  temperature  difference  has  increased. 
Fxsewhere,  very  little  has  changed.  Small-scale  irregularities  have  been 
smoothed  slightly  from  the  effect  of  horizontal  eddy  diffusion  in  the 
model.  Slight  changes  from  advection  as  well  as  diffusion  can  be  seen 
in  the  Kuroshio  and  East  Australia  Current  regions.  The  temperature 

pattern  southwest  of  Australia  has  changed  as  a res  lit  of  the  unrealistic 
current  there. 

The  velocity  fields  at  the  end  of  day  105  are  shown  in  Fig.  12. 

Again,  most  o’  the  changes  have  occurred  in  low  latitudes,  and  perhaps 
the  most  striking  feature  is  the  zonal  flow  that  has  developed  near  the 
equator.  In  rough  agreement  with  nature,  the  flow  is  directed  westward 
in  the  upper  layer,  eastward  in  the  lower  (~  Cromwell  Current)  la'er. 
However,  the  speed  in  the  lower  layer,  10  cm/sec,  is  small  by  a con- 
siderable margin  compared  with  Cromwell  Current  speeds  of  order  100  cm/sec. 
The  discrepancy  arises  at  least  partly  because  the  velocities  are  layer- 
average  values  (cf.  Fig.  1).  Thus  1C  cm/sec  applies  to  the  entire  lower 
layer  of  thickness  200  m. 

The  flow  in  the  California  and  Peru  Current  regions  remains  confus- 
ing. In  addition,  there  are  no  North-  and  South-Equatorial  counter- 
currents  (eastward-directed)  in  the  upper  layer.  This  may  be  partly 
related  to  a relatively  large  horizontal  grid  spacing. 


Fifteen-day  prediction  of  temperature  fields  (end  of  day  105)  subject  to 
the  heating  of  Fig. 7b.  Contour  intervals  and  dashed  lines  same  as  Fig. 6. 
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ower  layer.  Note  the  more 
equator  at  both  levels. 
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VII.  DISCUSSION  AND  CONCLUS IONS 


The  results  from  the  Pacific  model  are  generally  encouraging.  Per- 
haps to  some  extent  they  are  also  surprising,  in  that  an  ocean  model 
can  apparently  make  15-dav  predictions  showing  appreciable  changes  (of 
order  1°C)  in  the  temperature  distributions,  at  least  in  equatorial 
regions.  The  thermal  inertia  of  the  ocean  is  large  compared  with  that 
of  the  atmosphere.  But  this  apparently  does  not  preclude  appreciable 
oceanic  changes  (in  the  model)  on  time  scales  of  order  10  days  for  the 
upper  layers  : n low  latitudes.  This  may  have  important  implications 

for  monthly  and  seasonal  simulations  with  interacting  atmospheric  and 
oceanic  models. 

There  are,  of  course,  a number  of  deficiencies  in  the  Pacific  model. 
For  the  most  part  these  might  have  been  anticipated  from  the  known  limi- 
tations of  the  model  design.  The  model  should  not  be  taken  too  seriously 
near  eastern  boundaries,  where  it  generally  predicts  doWUelling.  This 
seems  to  be  a characteristic  shortcoming  of  general  circulation  models 
with  coarse  grid  spacing  (e.g.  , Bryan  and  Cox  1967).  It  is  likely  that 
inclusion  of  continental  slope  and  shelf  topography,  as  well  as  much 
finer  horizontal  grid  resolution,  will  be  required  before  models  such 

as  this  can  adequately  simulate  upwelling  near  the  western  shores  of 
continents . 

Regarding  the  confusing  flow  in  the  eastern  tropical  region,  it  is 
less  clear  what  the  cause  might  be.  It  is  doubtful  that  the  southward- 
flowing ( El  Nino)  current  can  be  believed  to  extend  to  northern  Chile, 
because  such  an  event  has  apparently  never  been  observed.  The  fault 
may  lie  either  with  the  model,  or  with  the  input  data  (wind  stress, 
heating,  and  initial  temperatures),  or  both.  For  example,  neither  the 
coarse  grid-resolution  of  the  model  nor  that  of  the  input  data  can  ade- 
quately resolve  the  smaller- scale  temperature  structure  observed  near 
Peru.  (See,  e.g.,  Wyrtki  1969.)  This  problem  is  at  least  partly  re- 
lated to  the  upwelling  difficulty. 

A related  limitation  is  that  the  model  does  not  predict  mixed-layer 
thickness,  but  rather  assumes  it  to  be  uniformly  50  m (one-half  the 


upper-layer  thickness  of  Fig.  1).  The  limitation  is  not  quite  that 
severe  because  the  model  can  predict,  e.g.,  0^  = 0^  at  some  point,  in 
which  case  the  mixed  layer  extends  locally  to  300  m depth.  Neverthe- 
less, this  is  one  of  several  limitations  motivated  largely  by  a desire 
for  simplicity.  A variable  mixed-layer  thickness  would  add  still  another 
prediction  equation  to  the  model.  The  strategy  has  been  to  try  to  ex- 
tend the  simpler  model  as  far  as  possible,  recognizing  that  improvements 
might,  and  probably  should,  be  added  later. 

As  a consequence  of  the  assumed  constant  levels  of  Fig.  1,  the 
choices  made  for  layer  thicknesses  represented  trade-offs  among  many 
conflicting  requirements.  For  example,  a horizontal  bottom  associated 
with  the  base  of  the  main  thermocline  in  middle  latitudes  might  be 
placed  more  realistically  at  600  m to  800  m depth.  However,  it  seemed 
more  important  to  attempt  a better  model  for  the  vast  tropical  region. 
Here,  the  relatively  shallow  Cromwell  Current  dictated  a bottom  depth 
of  perhaps  300  m for  a two-level  model.  The  latter  was  the  decisive 
factor  in  this  instance.  However,  the  assumption  of  zero  vertical 
velocity-shear  in  the  lower  half  of  the  lower  layer  (cf.  Fig.  1)  is 
perhaps  more  of  a compromise  in  favor  of  middle  latitudes.  In  the 
Cromwell  Current,  appreciable  vertical  shears  extending  to  300  m depth 
and  beyond  have  been  reported  by  Knauss  (1960)  and  others. 

In  spite  of  these  and  similar  criticisms  of  the  model,  the  results 
on  balance  are  encouraging.  The  ocean  model  with  its  coarse  grid  may 
be  unable  to  resolve  the  small-scale  temperature  structure  off  Peru; 
but  neither  could  a (hypothetically)  coupled  atmospheric  general  circu- 
lation model  with  its  even  larger  grid  spacing.  Perhaps  the  most  en- 
couraging factor  arises  from  an  apparent  paradox  in  the  results:  we 

input  observed  data,  initialized  the  model,  and  then  made  a "forecast." 

If  the  data  and  the  model  were  both  "correct,"  then  the  forecast  might 
have  properly  been  for  no  change.  However,  the  real  ocean  is  probably 
never  in  equilibrium.  Moreover,  there  is  no  obvious  reason  to  expect 
consistency  among  an  ocean  model  or.  the  one  hand,  and  distributions  of 
temperature,  wind  stress,  and  heating  on  the  other,  all  produced  by 
different  investigators  working  independently  of  each  other.  Our  input 
temperature  data  were  apparently  incomplete  and  inconsistent  with  the 
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model  in  the  equatorial  region.  The  encouraging  factor  is  that  the 
model  predicted  the  cold  tongue  of  water  toward  the  east,  a prediction 
that  seems  to  be  in  better  agreement  with  observation,  and  which  is  of 
sufficient  extent  that  it  might  be  detectable  by  an  atmospheric  circu- 
lation model. 

None  of  the  foregoing  is  meant  to  suggest  that  ..ie  prediction  of 
a cold  geographic  equator  and  zonal  flow  is  anything  new.  Haney  (1971a) 
found  a similar  result  with  a filtered,  six-level  model  of  a rectangular 
ocean.  One  point  of  distinction  concerns  time  scales,  which  are  vastly 
different:  200  to  210  years  in  Haney's  case  (the  model  times  for  which 

Haney  presents  results)  versus  90  to  105  days  in  the  present,  unfiltered 
model.  The  purposes  of  the  two  models  are  also  different.  Haney's, 
much  like  the  early  models  of  Bryan  and  Cox,  has  the  major  purpose  of 
studying  the  long-term,  steady-state  response  of  the  oceans.  His  model 
is  more  satisfying  in  the  sense  that  it  is  able  to  predict  features 
starting  from  simpler  initial  conditions,  e.g. , isothermal  and  a state 
of  rest.  The  present  model  cannot  do  this  because  its  time  increments 
are  much  too  small.  Rather,  its  stated  purpose  is  to  make  short-term, 
monthly  and  seasonal  predictions,  and  the  present  model  must  be  initial- 
ized from  observed  data.  Thus,  in  another  sense  the  present  model  may 
be  more  satisfying  as  a shorter-term  prediction  model.  The  point  is 
that  the  predicted  rapid  response  of  order  10  days  in  low  latitudes 
appears  to  be  new. 

The  predicted  zonal  flow  and  temperature  anomaly  at  the  equator 
(anomalous  from  an  atmospheric  point  of  view)  might  be  further  criti- 
cized as  being  too  special.  Although  the  stated  purpose  of  this  work 
is  to  develop  a world  ocean  model  for  climate  experiments , and  not  to 
investigate  equatorial  effects  as  such,  there  are,  nevertheless,  valid 
questions.  How  sensitive  is  the  result  to  the  initial  temperature  dis- 
tribution? or  to  the  applied  wind  stress?  or  to  the  basin  geometry? 

The  foregoing  result  was  in  fact  found  to  be  typical  of  several 
runs  made  with  the  Pacific  model,  and  of  many  runs  made  with  earlier, 
simplier  versions.  The  simplest  case  was  that  of  a (two-level)  homo- 
geneous ocean  in  a rectangular  basin  subject  to  a simple  zonal  wind 
stress  distribution  having  a westward  component  at  the  equator.  Com- 
parison runs  fc  different  vertical  viscosities  and  layer  thicknesses 
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all  showed  the  same  qualitative  results  near  the  equator.  The  zonal 
components  of  flow,  westward  (eastward)  at  the  upper  (lower)  level, 
were  often  realistic.  However,  the  meridional  components  were  gen- 
erally at  least  as  large,  and  these  were  poleward  (equatorward)  at  the 
same  respective  levels.  (Pure  zonal  flow  can  of  course  occur  at  the 
equator  for  north/south  symmetry  of  tha  boundary  conditions.) 

Qualitatively  similar  results  were  obtained  in  the  simple  baro- 
clinic  comparison  runs — but  only  for  the  initialization  state--whenever 
the  initial  temperature  was  a specified  function  of  latitude  and  depth, 
but  not  of  longitude.  Predictions  then  showed  the  simultaneous  develop- 
ment of  more  realistic  zonal  flow  near  the  equator  together  with  an 
east-west  temperature  gradient,  with  colder  water  toward  the  east. 
However,  well-defined  countercurrents  displaced  north  or  south  of  the 
equator  were  generally  not  found  in  the  results.  This  difficulty  was 
also  encountered  by  Haney  (1971a)  with  a six-level  model.  The  point 
here  is  that  the  foregoing  results  in  the  equatorial  zone  were  not  an 
isolated,  special  case.  The  results  suggest  that  the  observed  zonal 
flow  near  the  equator  and  the  temperature  anomaly  toward  the  east  may 
be  intimately  coupled.  However,  a physical  explanation  is  not  offered 
here,  and  further  study  with  a more  specialized  equatorial  model  is 
required  to  clarify  this  point. 

A f nal  encouraging  factor  from  the  Pacific  model  is  the  speed  of 
the  calculations.  As  mentioned  at  the  end  of  Sec.  V,  approximately 
4 min  of  IBM  360/91  machine  time  are  required  per  simulated  day  of  the 
Pacific  model  on  a 2 x 2 deg  mesh.  Doubling  this  gives  approximately 
8 min/day  for  the  global  case,  a conservative  estimate  providing  high- 
latitude  boundaries  are  maintained  for  linear  computational  stability. 
Increasing  the  2x2  deg  mesh  size  to  the  4x5  deg  global  mesh  in  cur- 
rent use  with  the  Mintz-Arakawa  model  gives  a coarse-grid  case  for 
direct  comparison.  Computer  times  in  the  coarse-grid  case,  taking  into 
account  the  allowable  increase  in  time  increment,  should  decrease  by 
a factor  of  ten  or  so.  This  gives  less  than  1 min/day  compared  with 
about  14  min/day  for  the  faster  version  of  the  Mintz-Arakawa  model. 

Thus,  computer  times  for  the  global  ocean  model  may  be  only  a tiny 
fraction  of  the  total  time  required  in  projected  atmosphere-ocean  model 
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experiments,  but  in  any  event  no  more  than  perhaps  one-third,  depend- 
ing on  the  mesh  size. 

In  conclusion,  the  Pacific  model  is  able  to  simulate  reasonably 
well  most  of  the  gross  features  of  the  circulation  and  temperature  dis- 
tribution in  the  upper  ocean.  The  model  is  able  to  predict  temperature 
changes  of  order  1°C  on  time  scales  of  order  10  days  in  low  latitudes, 
indicating  that  projected  monthly  and  seasonal  simulations  with  a coupled 
atmosphere-ocean  model  may  be  productive.  Finally,  computer  costs  appear 
modest  for  an  unfiltered  model.  The  Pacific  model  is  now  being  extended 
to  the  world  ocean  case,  with  possible  improvements  postponed  to  later. 
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